Minimum energy paths for dislocation nucleation in strained epitaxial layers 
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We study numerically the minimum energy path and energy barriers for dislocation nucleation in 
a two-dimensional atomistic model of strained epitaxial layers on a substrate with lattice misfit. 
Stress relaxation processes from coherent to incoherent states for different transition paths are 
determined using saddle point search based on a combination of repulsive potential minimization 
and the Nudged Elastic Band method. The minimum energy barrier leading to a final state with 
a single misfit dislocation nucleation is determined. A strong tensile-compressive asymmetry is 
observed. This asymmetry can be understood in terms of the qualitatively different transition paths 
for the tensile and compressive strains. 
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The growth and stability of heteroepitaxial overlay- 
ers is one of the central problems in current materials 
science. Energy-balance arguments for the competition 
between strain energy build-up and strain relief due to 
dislocation nucleation in mismatched epitaxial films lead 
to the concept of an equilibrium critical thickness. This 
is defined as the thickness at which the energy of the epi- 
taxial state is equal ±o that of a state containing a sin- 
gle misfit dislocationEl. The predicted critical value from 
this equilibrium consideration however, both from con- 
tinuous elastic madelsatl and from models incorporating 
layer discretenessH, is much smaller than the observed ex- 
perimental value for the breakdown of the epitaxial state. 
This suggests that the defect-free (coherent) state above 
the equilibrium critical thickness is metastableu and the 
rate of dislocation generation is actually controlled by ki- 
netic considerations. The idea of strain relaxation as an 
activated process is supported by experimental resultsii 
the temperature dependence of the critical thicknessli 
ft is also the fundamental assumption in kinetic semi- 
empirical modelaJ. 

Physically, it is expected that the lowest energy barrier 
for the nucleation of dislocations would correspond to a 
path that initiates from the free surface (with or with- 
out defects). Such processes have been con^dered in a 
number of studies within continuum modelsutj. It has 
been pointed out that surface steps and surface rough- 
ness that are not included in the continuum model iCOiiJd 
play an important role for dislocation nucleationOila. 
Thus, atomistic study is important for a detailed under- 
standing and direct determination of the mechanisms for 
defect nucleation in epitaxial films. However, determina- 
tion of the correct transition path and the nucleation bar- 



rier from the initial coherent state to the final state with 
misfit dislocations is an extremely challenging problem in 
an atomistic model. There are many saddle points and 
transition paths in the neighborhood of the initial coher- 
ent state. A brute force molecular dynamics (MD) study 
is impractical because of the rare event nature of the acti- 
vated processes. In recent years, great progress has been 
made in the general field of sea rch f or transition paths for 
complicated energy landscapestalla j-Jp^addition, various 
accelerated hyperdynamics schemesEJ'EEl have been devel- 
oped to overcome the computational problems for rare 
events. However, these schemes still involve consider- 
able computational efforts for complicated, large energy 
barriers and often require a qualitative picture of the en- 
ergy landscape as a starting point. Recently we hawe 
developed a repulsive potential minimization methodtj 
which allows automatic generation of many paths lead- 
ing away from the initial minimum energy coherent state. 
When .combined with the Nudged Elastic Band method 
(NEB)lia, this method can be used to efficiently locate 
saddle point configurations and barriers for each transi- 
tion path without having to make any specific assump- 
tions about the nature of the transition path. 

For the present study, we consider a two-dimensional 
model of the epitaxial film and substrate where the 
atomic layers are confined to a plane as illustrated in Fig. 
|l| (a). Interactions between atoms in the system are mod- 
elled by a generalized Lennard-Jones pair potentialElj, 
that is modified ta insure that the potential and its first 
derivative vanisbO at a cut-off distance r c as 

U(r) = V(r), r<r ; 
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U(r) = V(r) 
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where 



V(r) = e 



r>r , (1) 



(2) 



and r is the interatomic distance, e the dissociation en- 
ergy and ro the equilibrium distance betaieen the atoms. 
This potential has been used previously^, with n = 12 
and m = 6, in a Monte Carlo simulation of epitaxial 
growth. We have chosen the value n = 8 and m = 5 
for the present study. In contrast to the standard 6—12 
potential, this 5 — 8 potential has a slower fall-off. Thus, 
when combined with the variation of the cutoff radius 
r c , this choice allows us to systematically study the ef- 
fect of the range of the potential on misfit dislocation. 
Also, the 5 — 8 potential gives a more realistic descrip- 
tion of metallic systems than the 6 — 12 case. The equilib- 
rium interatomic distance ro was set to different values 
for the substrate, epitaxial film, and the substrate-film 
interfaces. The substrate ro = r ss and the epitaxial film 
To = r ff parameters were varied to give a misfit / be- 
tween lattice parameters as / = (rg — r ss )/r ss . For the 
film-substrate interaction we set the equilibrium distance 
as the average of the film and substrate lattice constants, 
Tq = f{ s = (rs + r ss )/2. Positive misfit / corresponds to 
compressive strain and negative / to tensile strain when 
the film is coherent with the substrate. Calculations were 
performed with periodic boundary conditions in the di- 
rection parallel to the film-substrate interface. Typically, 
one-dimensional layers containing 50 atoms or more were 
used in the calculations. In the calculations the bot- 
tom five layers represented the substrate, with the last 
two layers held fixed to simulate a semi-infinite substrate 
while all other layers were free to move. 

Our new scheme of identifying the saddle points and 
the transition paths consists of several stages. First, the 
initial epitaxial state is prepared by minimizing the to- 
tal energy of the system using MD cooling. This leads 
to an initial coherent epitaxial state in which the inter- 
layer spacing is relaxed, but the atoms within the layers 
are under compressive or tensile strain according to the 
misfit. Next, we introduce a short-ranged repulsive po- 
tential centered at the starting epitaxial configuration of 
the form 



Utot{r) = U(r) + Aexp{-a(r - r ) 2 }, 



(3) 



where ro's a&e, the coordinates of the initial state at 
the minimumta. The basic idea here is to modify the 
local energy surface to make the initial epitaxial state 
unstable. When the system is slightly displaced from 
the initial state (randomly or in a selective way), it 
will then be forced to move to nearby minimum energy 
states. By making this repulsive potential sufficiently 
localized around the initial potential minima, the sur- 
rounding minima would be unaffected and the final state 



energy would then depend only on the true potential of 
the system. This method can generate many different 
final states depending on both the initial displacements 
and the parameters of the repulsive potential. In this 
work, we only consider final configurations corresponding 
to a single misfit dislocation. Rather than trying random 
initial displacements, some knowledge of the dislocation 
generation mechanism is useful for expediting the pro- 
cess. Given the knowledge of the final state, we then use 
the NEB method to locate the saddle point and find the 
minimum energy path (MEP) between the initial and fi- 
nal states. We note that the path generated in the first 
part of moving away from the repulsive potential can be 
used as an initial guess for the MEP determination in the 
NEB method. 

We have performed calculations for various misfits but 
in this paper we concentrate on the case of a relatively 
large 8% misfit. We work with dimensionless quantities 
with e as the energy scale and r ss as the length scale. 
Two different choices of cutoff, namely r c =1.5 r ss and 
r c =4.7 r ss were made. The results for the barriers are 
qualitatively similar, so we present here only the results 
for the short range potential with r c =1.5 r ss . However, 
the equilibrium critical thickness and its asymmetry with 
respect to tensile or compressive strain are very sensitive 
to the range of the potentiaO. 

The results for the MEP from coherent to incoherent 
states are shown in Fig. [I] for a film under tensile strain 
and Fig. ^| for compressive strain. They show clearly 
the existence of an energy barrier for the nucleation of 
a misfit dislocation. Thus, the non-equilibrium critical 
thickness can be much larger than the equilibrium value 
and it is controlled in practice by the kinetics of defect 
nucleation. 

For compressive strain, the final state is characterized 
by the presence of an adatom island on the surface of the 
film for each misfit dislocation. The number of adatoms 
in the island exactly corresponds to the number of layers 
in the film. Such form of the final state is determined by 
the geometry of the misfit dislocation. For every misfit 
dislocation, an extra atom is removed from each layer to 
relieve the compressive stress. For tensile strain, the fi- 
nal state is characterized by the presence of pits on the 
surface. Again, the size of the pit is determined by geo- 
metrical considerations. For every misfit dislocation, an 
extra atom has to be added to each layer to relieve the 
tensile stress. For both cases, the dislocation core is lo- 
calized in the substrate-film interface region. 

Figures [l] and || also show the particle configurations at 
the different points along the MEP which reveal details 
of defect nucleation and strain relaxation process. The 
transition path for the compressive strain has a more lo- 
cal nature, with relatively fewer bonds involved initially, 
whereas for the tensile strained film, the nucleation pro- 
ceeds via a more collective path, involving concerted mo- 
tion along glide planes. The energy barrier for nucleation 
of a dislocation is much higher for the compressive strain 
relative to the case of tensile strain. This asymmetry is 
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very robust and it persists when we change the range of 
the potential by varying the cutoff. 

To understand the origin of this asymmetry, we plot 
in Fig. H the distribution of the nearest-neighbor bond 
lengths for the film from the initial epitaxial film to 
the saddle point configuration for both the compressive 
and tensile cases. It can be seen that the behavior of 
the compressively strained film and the tensile-strained 
film is very different. In the tensile case, the redistribu- 
tion of the bond lengths going from the initial coherent 
state to the saddle point configuration involves a signif- 
icant contraction of the intralayer bonds leading to par- 
tial relaxation of the tensile strain in the film. On the 
other hand, for the compressively strained film, the ini- 
tial delta function peak for the intralayer bond lengths 
broadens almost symmetrically and there are no signifi- 
cant relaxation of the compressive strain in the film. This 
explains the relatively higher energy costs and a corre- 
sponding larger nucleation barrier for the compressive 
strained film. Microscopically, the origin of the differ- 
ent behavior could arise from the strong anharmonicity 
of the interaction potential. For the compressive strain, 
intralayer rearrangements involve some further compres- 
sion of the bonds which is energetically costly. Thus, a 
more localized initial configuration with a higher barrier 
results as opposed to the collective behavior of the ten- 
sile strained layer. We have also checked that the bound- 
ary conditions and system sizes do not affect the results 
qualitatively by comparing results from systems with pe- 
riodic and free boundary conditions, and for layers twice 
as long. 

In summary, we have developed a new scheme of iden- 
tifying minimal energy path for spontaneous generation 
of misfit dislocation in an epitaxial film. This new ap- 
proach requires no a priori assumptions about the nature 
of the transition path or the final states. A nonzero ac- 
tivation barrier for dislocation nucleation is found in the 
minimum energy path from coherent to incoherent state 
above the equilibrium critical thickness, confirming the 
meta-stability of the epitaxial coherent film. The nucle- 
ation mechanism from a flat surface depends crucially 
on whether we start from a tensile or compressive initial 
state of the film. This asymmetry originates from the an- 
harmonicity of the interaction potentials which leads to 
qualitatively different transition paths for the two types 
of strains. A tensile-epmpressive asymmetry has also 
been found previously! 1 in other contexts. The present 
method can be extended to three-dimensional models 
with more realistic interaction potentials. Preliminary 
calculations for the Pd/Cu and Cu/Ed systemsEj with 
the Embedded Atom Model potentials^ confirms the ef- 
fectiveness of the method in three dimensions. These 
results will be published elsewhere. 
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FIG. 1. Particle configurations and energy change _E, — Eq 
at different states (images) along the minimum energy path, 
for tensile strain (/ = —8%). The two layers at the bottom 
are held fixed while all others are free to move. Open circles 
represent the substrate atoms and filled circles the epitaxial 
film. Only the central part of the layers with major atom 
rearrangements is shown. 

FIG. 2. Same as Fig. |l| but for compressive strain 
(/ = +8%). 
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FIG. 3. Nearest-neighbor bond distributions of the epitax- 
ial film at the saddle point for the (a) tensile, and (b) com- 
pressive cases. Solid and dotted arrows indicate the position 
of the delta-function peak corresponding to intralayer and in- 
terlayer bond distributions of the initial coherent film. 
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